

In [1]:

    
import numpy as np
import matplotlib.pyplot as pl
import seaborn as sn
from sklearn.cluster import MiniBatchKMeans
from sklearn import metrics
%matplotlib inline



In [37]:

    
stI=np.load('../data/stokesI_sunspot.npy')
stV=np.load('../data/stokesV_sunspot.npy')
wave = np.loadtxt('../data/wavelengthHinode.txt')
stI.shape









    Out[37]:





(512, 512, 112)

Index

Contrast and velocity fields

Classification

Contrast and velocity fields



In [3]:

    
sn.set_style("dark")
f, ax = pl.subplots(figsize=(9,9))
ax.imshow(stI[:,:,0], aspect='auto', cmap=pl.cm.gray)









    Out[3]:





<matplotlib.image.AxesImage at 0x7728990>






    



/usr/pkg/python/Python-2.7.6/lib/python2.7/site-packages/matplotlib-1.4.3-py2.7-linux-x86_64.egg/matplotlib/font_manager.py:1282: UserWarning: findfont: Font family [u'Arial'] not found. Falling back to Bitstream Vera Sans
  (prop.get_family(), self.defaultFamily[fontext]))

Let us compute simple things like the contrast and the Doppler velocity field



In [9]:

    
contrastFull = np.std(stI[:,:,0]) / np.mean(stI[:,:,0])
contrastQuiet = np.std(stI[400:,100:300,0]) / np.mean(stI[400:,100:300,0])
print("Contrast in the image : {0}%".format(contrastFull * 100.0))
print("Contrast in the quiet Sun : {0}%".format(contrastQuiet * 100.0))









    



Contrast in the image : 24.5026469231%
Contrast in the quiet Sun : 6.59198462963%

Now let us compute the velocity field. To this end, we compute the location of the core of the line in velocity units for each pixel.



In [31]:

    
v = np.zeros((512,512))
for i in range(512):
    for j in range(512):
        pos = np.argmin(stI[i,j,20:40]) + 20
        res = np.polyfit(wave[pos-2:pos+2], stI[i,j,pos-2:pos+2], 2)       
        w = -res[1] / (2.0 * res[0])
        v[i,j] = (w-6301.5) / 6301.5 * 3e5



In [34]:

    
f, ax = pl.subplots(figsize=(9,9))
ax.imshow(np.clip(v,-5,5))
f.savefig('velocities.png')



In [6]:

    
f, ax = pl.subplots(nrows=1, ncols=2, figsize=(15,9))
ax[0].imshow(stI[:,0,:], aspect='auto', cmap=pl.cm.gray)
ax[1].imshow(stV[:,0,:], aspect='auto', cmap=pl.cm.gray)









    Out[6]:





<matplotlib.image.AxesImage at 0x8379150>



In [8]:

    
f.savefig('exampleStokes.png')

Classification



In [17]:

    
X = stV[50:300,200:450,:].reshape((250*250,112))



In [22]:

    
maxV = np.max(np.abs(X), axis=1)
X = X / maxV[:,None]



In [40]:

    
nClusters = 9
km = MiniBatchKMeans(init='k-means++', n_clusters=nClusters, n_init=10, batch_size=500)



In [41]:

    
km.fit(X)









    Out[41]:





MiniBatchKMeans(batch_size=500, compute_labels=True, init='k-means++',
        init_size=None, max_iter=100, max_no_improvement=10, n_clusters=9,
        n_init=10, random_state=None, reassignment_ratio=0.01, tol=0.0,
        verbose=0)



In [42]:

    
out = km.predict(X)



In [43]:

    
avg = np.zeros((nClusters,112))
for i in range(nClusters):
    avg[i,:] = np.mean(X[out==i,:], axis=0)



In [44]:

    
f, ax = pl.subplots(ncols=3, nrows=3, figsize=(12,9))
loop = 0
for i in range(3):
    for j in range(3):
        percentage = X[out==i,:].shape[0] / (250*250.) * 100.0        
        ax[i,j].plot(km.cluster_centers_[loop,:])
        ax[i,j].set_title('Class {0} - {1}%'.format(loop, percentage))
        loop += 1
pl.tight_layout()



In [29]:









    Out[29]:





(9, 112)



In [ ]: